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Abstract. Macromolecular complexation leading to coupling of two or more cellular 
membranes is a crucial step in a number of biological functions of the cell. While 
other mechanisms may also play a role, adhesion always involves the fluctuations 
of deformable membranes, the diffusion of proteins and the molecular binding and 
unbinding. Because these stochastic processes couple over a multitude of time and 
length scales, theoretical modeling of membrane adhesion has been a major challenge. 
Here we present an effective Monte Carlo scheme within which the effects of the 
membrane are integrated into local rates for molecular recognition. The latter step 
in the Monte Carlo approach enables us to simulate the nucleation and growth of 
adhesion domains within a system of the size of a cell for tens of seconds without loss 
of accuracy, as shown by comparison to 10® times more expensive Langevin simulations. 
To perform this validation, the Langevin approach was augmented to simulate diffusion 
of proteins explicitly, together with reaction kinetics and membrane dynamics. We use 
the Monte Carlo scheme to gain deeper insight to the experimentally observed radial 
growth of micron sized adhesion domains, and connect the effective rate with which 
the domain is growing to the underlying microscopic events. We thus demonstrate that 
our technique yields detailed information about protein transport and complexation 
in membranes, which is a fundamental step toward understanding even more complex 
membrane interactions in the cellular context. 
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At the origin of many biological phenomena is cell adhesion promoted by the formation 
of macromolecular ensembles. Despite intensive research over the last two decades 
and the pressing biological signihcance [I4| - [I^ , the growth of these structures 
in membranes is still poorly understood. Formation of adhesions involves a number 
of stochastic events occurring on different length and timescales. The minimal system 
involves protein diffusion and formation of bonds, which occurs on characteristic times 
of 10“^ —10“^ s. These two processes couple to fast membrane fluctuations (10“® — 10“® 
s). Several length scales are also involved - from nanometer separations necessary for 
molecular recognition to the micron-sized objects that are being grown. Moreover, 
molecular complexation induces membrane deformations which in turn promotes long- 
range cooperative effects. If all these elements are considered, difficulties in modeling 
the dynamics of macro-molecular scaffolding come as no surprise. 

Early attempts to model the formation of macromolecular structures were related 
to interactions of protein-decorated membranes with underlying substrates containing 
the appropriate binding partners in the adhesion process. Thereby, analogies with 
classical theories of growth (Stefan problem and kinetically limited aggregation) were 
explored 19 22 . Other approaches focused on the role of the membrane fluctuations 


23 ,^. Furthermore, a number of scaling laws were suggested after the analysis of 
the relationship between the various involved stochastic processes 25 . However, only 


limited experimental conhrmation has been obtained to support these arguments 26 27 


Later efforts concentrated on the construction of accurate simulation schemes that 
treat the membrane fluctuations explicitly. First, dynamics of domain formation was 
studied by Monte Carlo approaches where furthermore the diffusion was treated by 
a random walk and complexation of proteins was explored through Metropolis rates 
j5,28,^. Concomitantly, Langevin simulations j4,30 33 were developed. In earlier 
attempts (4 |[^[m], binding and unbinding was not considered, while latter efforts 

The 


involved rates that are functions of the instantaneous membrane prohle l32, 33 


problem with all these methods is that only micron-sized systems could be studied 
for about a millisecond. Consequently, long time-scale dynamics associated with the 
formation of larger macromolecular structures, such as radially growing domains and 
diffusion-limited aggregation, remained out of reach. To address these biologically 
relevant issues, signihcant efforts went toward developing coarse-grained simulation 
methods. This resulted in mapping the problem onto lattice gas and Ising-like models 


34 38 , which is, however, accurate only in a limited range of parameters. 


Here we build on the experience in coarse-graining the dynamics of nucleation 
of macromolecular complexes in membranes (^. We solve the problem of coupling 
time and length scales by constructing an effective Monte Carlo simulation scheme, for 
which we demonstrate applicability in a very broad range of parameters. The stepping 
stone for our approach is the realization that there is a clear separation of time scales 
between membrane fluctuations and protein binding and diffusion. This allows us to fully 
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circumvent simulating the membrane, by incorporating its influence into effective rates 
for the (de)complexation of proteins. We validate our scheme against explicit Langevin 
simulations which themselves were shown to agree very well with experiments 


in the context of the nucleation 40 and the morphology of adhesion domains 32 


In order to make this comparison easier, we first present the underlying theoretical 
model, its direct implementation into the augmented Langevin scheme and, then, the 
upscaling and the construction of the effective Monte Carlo scheme. Furthermore, we 
demonstrate the potential of the Monte Carlo scheme by simulating radially growing 
domains containing up to 10® ligand-receptor bonds over several seconds, as observed 
in analogous experiments. This allows us to explore the membrane associated processes 
with very high precision, and to provide deeper understanding of the overall dynamics. 


2. Model 



Figure 1. In the model, the fluctuating membrane carries mobile ligands, which 
bind to immobile receptors placed equidistantly on the surface (characteristic spacing 
d). The formation of bonds is associated with the deformation of the receptor and 
the membrane, the latter being subject to a nonspecific harmonic potential with a 
minimum at /ig. 

Our model system (see figure consists of a flexible membrane that is positioned 
above a solid substrate. Receptors on the substrate can form bonds with ligands, 
embedded in the membrane. The receptors are placed on a regular square grid and 
are immobile in the current context. The ligands can diffuse within the membrane until 
a bond is formed and therefore the membrane is locally pulled towards the substrate. 
More elaborate versions of our model allow for both binders to be mobile and coupled 
to different reservoirs, simulating a finite vesicle, or an infinite bilayer. Furthermore, 
binder species with different properties can be simultaneously introduced. 
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The membrane is described as a thin sheet with an energy given by the Helfrich- 
Hamiltionian IdTl 


nuMr)] = I d^r (^^{Ah{r)f + | [h(r) - hof^ . 


( 1 ) 


Specifically, the lipid bilayer is parametrized in the Monge gauge, where the height h{r) 
is given as a function of the position r of the membrane above the substrate (Fig. [^. A 
list of the variables and parameters in equation ([^ can be found in table Specifically, 
the first term in equation ([^ is the deformation energy of the membrane, that is itself 
a product of the bending rigidity k and the local mean curvature of the membrane. 

While the specific protein molecules embedded in the cell wall (or membrane) are 
usually considered to be responsible for cell adhesion, over the past two decades a 
realization emerged that the cell membrane itself, being a floppy sheet, adds another 
unavoidable, yet not fully understood, interaction with the opposing surface it binds to. 
Although this interaction does not depend at all on any specific proteins, it can have a 
major impact on the protein-mediated adhesion and can be viewed as a mechanism 
that controls the binding affinity to the cell-adhesion molecules |^. Such steric 
interactions typically maintain the two membranes at relatively large separations ho, 
which can be modeled by introducing a nonspecific harmonic potential of a strength 7 
with the minimum at ho |33|[4^[45] . The strength of this potential depends directly 
on the average intensity of membrane fluctuations that are themselves regulated by the 
tension in a membrane but also by numerous other factors such as the thickness and the 
composition of the glycocalyx. In the mimetic systems, this contribution is dominated 
by continuous interactions between the membrane and the substrate, such as gravity, 
Helfrich-repulsion, or Van-der-Waals forces 44,46 . The strength of this potential can 


be obtained experimentally by the analysis of membrane fluctuations 147,48 . 


2.2. The bonds 

We assume that the receptors are thermalized springs with stiffness A and rest length Iq. 
This leads to the following expression for the energy of the N), bonds in the membrane 

Nb 


^B[hir)] = ^5(r-ri) 


i=l 


^(h(r) -lof-Cb 


( 2 ) 


Here, Cb accounts for the bond enthalpy gain for forming a bond and 5(r — ri) is the 
Dirac-Delta function for the positions ri of the bonds. 

As we assume that the structural fluctuations of free receptors occur on a faster 
time scale than the membrane dynamics, each bond fulfills a local detailed balance for 
the transitions between the bound and unbound states, given by the rates k°^{h{r,t)) 
and (h(r, t)) as 


^off t)) 

^on t)) 


= exp 




loY 


(-b 


-In 

2 


Xa^ 


( 3 ) 
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Table 1. Variables and parameters of our Helfrich-Hamiltionian (see equations Q 
and Q). 


Quantity 

Meaning 

Unit 

a 

lattice constant 

10 nm 

ksT 

thermal energy at 300 K 

4.14 X 10-21 J 

K 

bending rigidity 

ksT 

h{r) 

membrane profile 

a 

7 

curvature of the interaction potential 

ksTja^ 

ho 

minimum of the interaction potential 

a 

Nb{t) 

number of bonds 

- 

Ti 

position of bond i 

a 

A 

stiffness of the bond/receptor 

keT/a^ 

h 

rest length of the bond/receptor 

a 

eb 

binding enthalpy 

ksT 


Following this condition, each bond is locally in thermal equilibrium with the 
instantaneous membrane profile. Here, a is the range of the interaction potential of the 
ligand-receptor bond and for simplicity, we set jS = {kBT)~^ = 1. Equation ^ depends 
on the stretching energy of the bond (first term in the exponent), the binding affinity 
(second term) and an entropic contribution (last term) which describes the suppression 
of fluctuations if a receptor is bound to a ligand. This entropic contribution lowers the 
effective binding affinity |4^ . 

Inspired by [8,49, ^ we choose the rates for the creation of a bond as 


/c“(h(r,f)) = fco\/^exp 


~{{h{r,t)-k)-a}'^ 


( 4 ) 


where ko is the intrinsic reaction rate. From this local on-rate and the detailed balance 
condition the local, off-rate can be determined readily 


r®(h(r,t)) 


ko exp [-Cb] exp 


A(h(r,t) 


Iq) a 


\a^ 

~Y~ 


( 5 ) 


which is proportional to the rate of the Bell-model j^, accounting for the force acting 
on a bond A (h(r, t) — Iq). 

We show the reaction rates eqs. (|^ and (|^ in Fig. The association rate 
adopts the form of a Gaussian with a width inversely proportional to the stiffness of the 
receptor. The maximal association rate is obtained at /q -I- a, which is at the outer edge 
of the potential well associated with an unperturbed receptor. The off-rate increases 
exponentially with the distance between the receptor and the ligand. Interestingly, if 
the ligand is in the middle of the binding region (Zq + «/2 ), the stiffness of the receptor 
does not affect the breaking of the bond (the bond is not stressed). 
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Figure 2. Local reaction binding (left) and unbinding (right) rates equations Q and 
([^ shown as function of the membrane height. Other parameters were set as in table 
except for the binding energy (et = 0). 


2.3. Diffusion 


Due to the membrane fluidity, the molecules within the bilayer diffuse on its surface j^. 
Even though there may be an influence of the membrane elasticity on the diffusion of 
embedded proteins (for example by the curvature that a protein induces in the membrane 
^), these effects seem to be small for experimental relevant parameters [54] . 
Therefore, we simulate the mobility of binders by a random walk, whereby two proteins 
interact laterally by a hard-core potential. The time step of the random walk is given 
by 


32 53 



( 6 ) 


with the diffusion constant D. In the current work, only the ligands embedded in the 
membrane of the vesicle are allowed to diffuse. However, it is straightforward to extend 
the simulation scheme to situations in which both binders retain lateral mobility and 
explore the surface of the membrane. The latter may be finite as in the case of vesicles 
and cells. These situations are simulated using periodic boundary conditions on the 
level of the system, with a selected area in the center of the simulation box representing 
the area of contact between two cells or the cell/vesicle and the substrate. Consequently, 
the formation of bonds can take place only within this region, and the remainder of the 
system will be depleted from the binders due to the accumulation in the zone of contact. 

Binders can be also embedded in bilayers, which provides a constant chemical 
potential. For simulations of interactions with vesicles, a contact zone is defined, and 
the periodic boundary conditions are imposed for the bilayer grid. However, to maintain 
the constant chemical potential (constant concentration of binders in the bulk), entering 
and exiting of a binder from the contact zone is associated with placing or removing a 
binder from a random position outside the contact zone. 
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Figure 3. Snapshot of a simulation run. The grey receptor can form bonds with the 
orange ligands embedded in the fluctuating membrane (blue). The ligands can diffuse 
freely within in the membrane, whereas the receptors are immobilized and placed on 
a square grid. 


3. Langevin simulation scheme 

3.1. Equation of motion for the membrane 

In this scheme, the membrane shape (see Fig. is determined explicitly in every time 
step. Thereby, the system is propagated in time by means of the Langevin equation in 
Fourier space (see e.g. j^|^) derived from the equations ([^ and (|^ 

= - A(k) I [nk^ + 7 ] (h(k, t) - 4,oAho) 



Here, A(k) is the Oseen tensor, describing the hydrodynamic interaction between 
membrane and surrounding fluid and A is the area of the membrane. The stochastic 
force .^(k) in the Langevin equation above is set by the temperature of the surrounding 
fluid. Thereby, the Oseen tensor is connected to the stochastic force by the fluctuation- 
dissipation theorem 

(e(k)e(k')) = 2fcsTA(k)5(k + k'). (8) 

The dehnition of the Fourier transformation of the Langevin equation is given by 

h(k) = I d^r exp (—ik ■ r)h(r); h(r) = ^ ^ exp (ik • r)h(k). (9) 

Ja a ^ 

In general, it was shown that the Oseen tensor depends on the geometry of the membrane 
[ 5 ^. However, this dependence is very weak for membranes far away from the substrate 
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and only relevant for the four largest modes of the membrane for the parameters used 
in the simulations. Thus, we use the Oseen tensor for a free membrane 


A(k) 


1 

Arjk' 


( 10 ) 


where rj is the viscosity of the surrounding fluid. The Oseen tensor for the k = 0 mode 
diverges. Following the Oseen tensor for this mode is set to 


A(k) 


Stit] 


( 11 ) 


The Langevin-equation Q is solved numerically with the help of the Euler-Maruyama 
scheme (see for example [^). The time step in this scheme has to be set below the 
smallest time scale of the membrane 


T 


{kn 


«^Lx + 7’ 


( 12 ) 


which is the typical relaxation time of the mode with the largest k in the simulation 
(fc ma x = y/^Ti) and is on the order of 10“® s. 


3.2. Simulation scheme 

The simulation is performed following the algorithm shown on the left in hgure 
The hrst step initializes the system. This involves the thermal equilibration of a free 
membrane obtained by executing 10® steps in the time loop explained below without 
the reactions and binder diffusion. After that, the ligands are placed randomly on their 
lattice, and the receptors are put on a grid of the second lattice. 

The second step is the initialization of the time loop, where the step accounts for the 
shortest characteristic membrane time scale {At = r(/cinax))- Every time step involves (i) 
the calculation of the force on the membrane induced by the formed bonds in real space; 
(ii) the transformation of this force to the Fourier space; and (iii) the determination of 
bending and unspecihc forces in Fourier space (first term in equation ([^). The sum of 
this forces is input to the Euler-Maruyama step, within which the membrane prohle is 
updated in Fourier space and transformed back to real space. This back-transformation 
is a prerequisite for the execution of the association and dissociation step. Here, the 
binding probabilities are obtained from the equations (|^ and ([^, in which the height of 
the membrane and the time step of the simulation are required. As the time scale of the 
reactions is much larger than the typical time scale of the membrane, these probabilities 
are rather small. 

Finally, the mobile ligands need to be displaced to one of the neighboring 
unoccupied sites. In principle, the diffusion of binders is characterized by the time step 
given by equation (§, in which case the probability to jump in any direction would be 
1/4. However, as the time scale of the diffusion is typically several orders of magnitude 
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Table 2. Parameters used in the simulations 


Parameter 

value 

a 

10 nm 

K 

10 ksT 

7 

3.125 X 10-3 keT/a^ 

A 

7.5 X 10-3 to 5 X lO-^fcsT/nm^ 

7] 

2.4 X lO-^fcBTsa-3 

ho 

8 a 

lo 

4 a 

d 

8 a 

Pi 

l/64a-3 

ko 

2 X lOMo 5 X 10V3 

db 

5 to lO/cfiT 

D 

5 X 10^ to 5 X lO^nm^/s 


larger than the step of the time loop (i. e. td S> ^(/Cmax)), fhe probability of a jump is 
rescaled to 


p = 




Atd 


(13) 


This new probability guarantees the correct diffusive behaviour of the ligands. 

In the current simulations, the ligands are immobilized after they form a bond with 
a receptor, which means that only free ligands diffuse. This restriction is motivated by 
the experimental observation that the bonds change position only if they are subject to 
a significant lateral force j^. After the diffusion has been resolved, a new iteration in 
the time loop is started, or the simulation is terminated. 

Computationally, most time in this simulation scheme is consumed by Fast Fourier 
Transformations of the membrane profile and the forces, which scale like A^log(A^), 
[N is the number of considered lattice points), and not linearly like other operations 
(diffusion and reaction kinetics). Furthermore, the time step has to be chosen very 
small to accurately describe the time evolution of the membrane, and a large number of 
replicas must be produced to obtain a statistically sound representation of the system. 
These are the main reasons which make this simulation scheme computationally very 
expensive allowing only for length scales of up to 1 /rm^ to be simulated for up to 10“^ s. 


4. Effective Monte Carlo simulation 

4-1- Effective rates 

The difficulties that arise with Langevin simulations could be circumvented if the explicit 
treatment of the membrane could be avoided. We achieve this goal in an effective Monte 
Carlo scheme which is based on the recently acquired understanding of the effects of 
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Binder diffusion 
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Figure 4. Simulation schemes. Left side: Langevin scheme. The membrane is 
simulated explicitly. However, computational expensive Fourier transformations have 
to be performed. Right side: Effective Monte Carlo scheme. The binding kinetics 
is simulated with the effective reaction rates (18 \ and circumventing the explicit 
treatment of the membrane. 
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the membrane on the formation of bonds 139,46,48 . This scheme relies on the fact that 


the typical time scale of the membrane fluctuations depends on the viscosity r] of the 
surrounding fluid 


^mpm 


4r]qo 2r]qo 


~ 2 X 10"^ 


(14) 


nq^ + j 7 

Here, go = ( 7 /^)^^^ is the inverse lateral correlation length for a membrane without 
tension. 

Importantly, even the slowest modes are signihcantly faster than the reaction 
kinetics for ligand-receptor binding (the fastest avidin-biotin in membranes was reported 

(§|42 


to take place at 


10^ s-i 


40,42)), while other pairs are found at 


102.-1 


Consequently, the membrane fluctuations can be regarded as equilibrated with hxed 
mean shape as long as the conhguration of bonds interacting with the membrane remains 
unchanged. During this time, the fluctuating membrane, and with it the ligands, sample 
the entire probability distribution of distances between ligands and receptors. In the 
following, we denote the height distribution at the considered binding site r before the 
bond has formed by p(h'’), and the height distribution after a bound ligand-receptor pair 
is formed by p{h^). The latter is non-trivial if the receptor or the bond itself maintains 
some flexibility. 

The hrst and the second moment of these typically Gaussian distributions can be 
calculated explicitly for an arbitrary bond conhguration |^. Specihcally, we calculate 
a functional integral over all membrane prohles weighted by their Boltzmann factor (see 


Appendix A for technical details). As result, we obtain the mean height 


{A(r)> ^ rpT) = 




J Eh kei (90 k - ril)GH(>')“‘ kei (go k - rjl) 


( 16 ) 


and the huctuation amplitude 

(/7r))-(/!(r))"^(g’'7'W = 

\ -1 

16kei(go |r - ri|) Gij(r)“^ kei (go |r-rj|) 
u 




( 16 ) 


TT^ 


The sum runs over all pairs of bonds in the membrane at the positions ri and rj, while 
kei(a;) is the Kelvin function 59 . The elements of the coupling matrix Gij{r) are the 


effects of the existing bonds on the shape and huctuations at the arbitrary position r, 
whereby the membrane mediated interaction between the bonds are comprised in the 


off-diagonal elements (see A.7 for the explicit form of the matrix). 

The average shape and the fluctuation amplitude of a membrane containing a small 
cluster of bonds are shown in the top panels of figure At large distances from the 
cluster, the membrane is on average flat since it resides and fluctuates in the minimum 
of the nonspecihc potential. Because of a relatively high concentration of bonds within 
the cluster the membrane is likewise flat on average, but much closer to the substrate. 
At the same time, its fluctuations are strongly suppressed. However, the shape and 
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Figure 5. Membrane profile (top, left), fluctuations (top, right), effective off-rates 
(bottom, left) and effective on-rates (bottom, right) for a given bond configuration 
(white bonds and black free binding sites). The bond stiffness A is set to infinity and 
the binding affinity £(, to zero for simplicity. The remaining parameters can be found 
in table [D 


fluctuations of the membrane are significantly different in the vicinity of the bonds at 
the edge and in the center of the cluster. 

We use the two height distribution functions to average the Bell-Dembo rates (eq. 
(|^ and ([^) at the position of a free or a bound receptor giving rise to effective binding 
and unbinding rates 


= j dhXh")r"(h'’) 


(17) 


Appropriately inserting equations 0 . 0 . ([I^, and (p!6| into the the above expression, 
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Ko^ 


yx 




\/27r(l + A (cT^)^) 


exp 


A [h^ - {a+ k)Y 

2(1 +A (a")") 



Iq) + a[A (cr^) 


1 ]}- + 


(18) 


In this manner, the effective rates describing the association and the dissociation of a 
bond depends on the exact position of a bond, and the time dependent configuration 
of all bonds in the system. Examples of such rates for one bond configuration can 
be seen in bottom panels of figure Obviously, the rates reflect the average shape 
and fluctuations within the membrane (top panels), which are the result of the bond 
configuration around the respective binding site. The dissociation rate of a bond at the 
rim is up to two orders of magnitude larger than for a bond deep within the domain (see 
figure]^. This is due to the stabilization effects of the neighboring bonds, which share 
the deformation load and cooperatively suppress the fluctuations. On the other hand, 
the association is the largest near the bond domain and exponentially decreases on the 
length scale of the lateral correlation length with increasing distance to the domain. 


4-2. Simulation scheme 

We can now construct a Monte Carlo simulation of the adhesion process in which only 
the reaction kinetics and the diffusion of binders must be treated explicitly. Thereby, 
the effective rates for breaking or forming a bond at the given site must be determined 
for each site in every time step. This in turns requires inverting the coupling matrix 
containing all bonds, for every site in every step. In order to make the simulation 
fast, we assume that only the first two sets of neighboring bonds affect the rates on 
a particular (un)binding site (figure]^. Consequently, only the configuration of bonds 
in the immediate environment is taken into account in the calculation of the effective 
rates. Since this environment consists only of 9 sites, all possible configurations can 
be explored a-priori, and their respective rates used to create a lookup table. This 
restriction to the next-nearest neighbours is justified because the binding rates decay 
very fast with increasing distance between the bonds. 

A flow chart of the Monte Carlo scheme is shown on the right panel of figure 
To initialize the system, all ligands and receptors are positioned on their respective 
grids as in the Langevin simulations (random and ordered distributions are possible). 
Furthermore, the characteristic time steps are determined. The time step of the 
simulation is given by the characteristic diffusion time At/) (equation (|^). The time 
step for the reaction kinetics is set to be = Ato/n, where n is the smallest integer 
satisfying the inequality < 1. From here the probabilities for binding and 

unbinding are calculated as K°^^°^AtB, and stored in a lookup table. 

The simulation step starts with the reaction loop which consists of n iterations. 
In each iteration, for every binding site (i) the bond configuration is determined, (ii) 
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Figure 6. Example for a possible bond configuration two (left) or three (right) 
bonds (red squares). The association (left) or dissociation (right) rate at the considered 
binding site (center square) is determined by identifying the bond configurations (red) 
around the binding site and retrieving the appropriate reaction rate from the lookup 
table. This is done for all binding sites during one iteration. 


the appropriate rate is retrieved from the lookup table, (hi) association or dissociation 
is attempted, and (iv) the bond conhguration is updated. Following the reaction loop, 
each binder attempts to move to a neighboring site in a same manner as in the Langevin 
scheme. This completes the simulation step and the system is propagated in time until 
the program is terminated. While the program allows for the diffusion of both binder 
types, the following discussion will be restricted to the case when the receptors are 
immobilized. 

The advantage of the Monte Carlo scheme is that it allows for a larger time step 
and avoids Fast Fourier Transformations limiting the Langevin code. This allows us to 
simulate length scales of several tens of micrometers and time scales of several seconds 
with the resolution of about 100 nm and 10'^ s, which is necessary to understand 
biological processes. 

5. Validation of the Monte Carlo scheme 

In order to evaluate the applicability of the effective rates, we perform an extensive 
comparison of the results of the Langevin and Monte Carlo simulations. For this 
purpose, all parameters, the system size, and the statistics of data acquisition in the 
two approaches is identical and no £t parameters are used in the following discussion. 

We explore a very wide range of parameters: from soft to rather stiff receptors, 
binding affinities from the unstable to the enthalpy dominated adhesion, fast and slow 
diffusion of ligands (equivalent to changing the attempt reaction frequency). 
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Figure 7. Distribution of nucleation times. The effective scheme, the Langevin scheme 
and the analytical model produce the same distribution of nucleation times (without 
fitting parameter). The analytical curve is determined from the equation (6) in Bihr 
et al. |39| . The intrinsic binding affinity for protein binding is set to ej, = 6.56 fc^T, 
while the diffusion constant is D = bpiM/s. All simulations were performed in a 
simulation box of 640 nm x 640 nm with the densities of receptors and ligands of 
Pr = Pb = 1.5625 X 10“^ nm“^. The intrinsic binding rate was set to kg = 10® s“®, and 
the receptors are modeled as springs of stiffness A = 2 x 10“^ ksTjwad'. 


Soft receptors Medium receptors Stiff receptors 



Figure 8. Time evolution of the number of bonds for A = 7.5 x lO^^fc^T/nm^ 
and fco = 1.6 x 10® s“®(left), A = 2 x lO^^fc^T/nm^ and fcg = 10® s“^ (middle) and 
A = 5 X kgT/xiwd and fcg = 6.1 x 10^ s“^ (right) (remaining parameters see table 
1^. We compare the effective scheme (full lines) with the Langevin scheme (dotted 
lines). 
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5.1. Early stages of domain formation - nucleation dynamics 

We first focus on the simulation of rare events such as is the nucleation of adhesion 
domains. The number of bonds in such a domain can be calculated explicitly within 
the capillary approximation Once this number is estimated we perform about 2000 
simulations with each method to generate the distribution of nucleation times (Fig. [^. 
Specifically, each simulation is set to start from an equilibrated box with zero bonds. 
When a cluster of bonds of critical size is formed anywhere in the system, the simulation 
is interrupted, and the time necessary to achieve this domain size is recorded. As shown 
in Fig. excellent correspondence of the coarse-grained and the higher-level simulation 
approach is obtained for the entire distribution of nucleation times. This agreement 
could have been anticipated from the successful comparison of Langevin simulations with 
the analytic model for the nucleation dynamics of a single seed, based on a simplified 
version of the here used effective rates |^. The current, more accurate approach fully 
validates the concept of the effective rates and enables studies of the early stages of 
the adhesion process in the regimes that are either not accessible to analytic modeling 
or are extremely demanding from the computational point of view. Examples of such 
regimes, which can be now addressed with ease, are fast nucleation, competitive growth 
of multiple seeds, or diffusion limited nucleation. 

5.2. Full dynamics 

Encouraged by our success reproducing the nucleation dynamics, we validate the 
Monte Carlo scheme by reproducing the results of the higher level scheme for the full 
dynamic adhesion process, i.e. nucleation, growth and saturation to equilibrium. More 
specifically, for each set of parameters we perform 200 runs over which we average the 
dynamic process. This level of accuracy was found previously to produce converged 
results for the Langevin scheme in thermal equilibrium }32|[M) . 

We first explore the correspondence of the two schemes when the diffusion of ligands 
is fast {D = 5 X lO’^nm^/s), for soft, moderately stiff, and stiff receptors (Fig. [^. In 
each graph, the number of bonds as a function of time is presented for three different 
binding affinities (the smallest being at the phase transition to the unstable adhesion 
dominated by unbinding, two intermediate affinities, and one large affinity where the 
unbinding is negligible). We find that except for critical fluctuations at the phase 
boundary (e^ = 6.79 (^ , the two approaches show extremely similar dynamics. 

Very similar results are obtained for slow diffusion of ligands (Fig. [^, where 
the adhesion dynamics is shown for three receptor rigidities, at an intermediate 
binding affinity. Equally good, quantitative agreement is obtained for all affinities 
above the transition energy (data not shown). These exceptional results validate the 
concept of effective rates, and establishes the Monte Carlo approach as a reliable and 
versatile method for the simulation of protein mediated membrane interactions. It 
should be noted that different effective Monte Carlo schemes, based on the integration 
of membrane fluctuations in the Hamiltonian were successful in comparison with 
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Figure 9. Simulation curves (full lines for effective scheme, dotted lines for Langevin 
scheme) for different values of x = A/(8 y^k 7) = 7.26fcBT and D = 5x 10®nm^/s), 

remaining parameters as in figure 


the Langevin simnlation [^, with the time scale of reactions being a free htting 
parameter. However, the accnracy of that scheme relied on the magnitude of the effective 
cooperativity parameter y; to be much smaller than one. This dimensionless parameter 
evaluates the fluctuations of the unbound membrane with respect to the fluctuations of 
free receptors 

A 

^ 

The accuracy of the current scheme does not depend on the effective cooperativity 
parameter, which for the systems shown in Fig 9 range from 0.53, for the softest 
receptors, to 3.54, for the stiff receptors. Actually, the regime of large effective 
cooperativity parameters seem to be very important in the context of experiments with 
cells or vesicles (^ . 



6. Simulations of radially growing domains 


One of the basic mechanisms for the growth of adhesion domains is their radial expansion 
from a stable nucleus. As observed both in the cellular and cell-mimetic context, 
with different ligand-receptor pairs, such growth occurs naturally in membranes where 
the characteristic nucleation time is small compared to the dynamics of the domain 
expansion, and is common in situations where one of the binding partners is immobilized 
58 ,^. Particularly well-studied are radially growing domains in ligand-decorated 
vesicles binding on a substrate functionalized with receptors 19,42 . In these systems. 


radial growth was used for the determination of the effective binding rate of various 
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ligand-receptor pairs. This rate was found to depend significantly on the properties of 
the membrane due to strong correlations between the bonds |4^ . 

The analysis of the growth dynamics reveals that the growth of the 

domain is diffusion limited and the area of the domain increases linearly in time if 
the concentrations of ligands is smaller than the concentration of receptors 21,42 


Otherwise, the growth is reaction limited, and the area grows quadratically. Treating 
the growth dynamics as a diffusion-reaction problem, the diffusion constant of ligands, 
and the effective binding constant was extracted from the data (^. However, very 
little is known about the relation of such macroscopic measurements with the underlying 
microscopic binding and unbinding events, as well as protein motions in the membrane. 

Unfortunately, the limited size of systems that can be studied with the Langevin 
scheme makes this approach unsuited for the analysis of the radial growth process. 
Nevertheless, using a large number of replicas to reconstruct the representative 
dynamics, effective affinity, as well as the growth patterns could be identihed in the 
reaction limited case 32 . However, the issue of the system size is particularly acute 


for the diffusion limited processes, when a depletion zone around the growing domain 
forms, and extends faster than the domain itself 19 21,42 . This regime, as well as 


the continuous dynamics in the reaction limited case can only be obtained with the 
effective Monte Carlo approach developed here. As it will be shown in this section, such 
a study should clarify how the cooperative effects transmitted by the membrane affect 
the microscopic rates and the overall dynamics. 


6.1. Simulation details 

We perform a series of Monte Carlo simulations, where we use two opposing square 
grids of a size of 40.96x40.96 pm^ in the diffusion limited case and of 10.24x 10.24 pm^ 
in the reaction limited case (typical sizes of a giant unilamellar vesicle). The hrst grid 
carries 2.5 x 10^ receptors (soft or stiff), immobilized on a lattice. To simulate diffusion 
or reaction limited growth, the second grid is decorated by randomly placing 5 x 10"^ 
diffusing ligands or placing immobile ligands above the receptors, respectively. These 
concentrations, as well as the other parameters parameters are strongly inspired by 
the analogous experimental realizations of the system 42 . Specihcally, the height 


of the membrane (hg — Iq = 55 nm), curvature of the nonspecihc potential (7 = 
3.125 X 10“^ ksT /a^), bending rigidity of the membrane {k = 10 ksT), binding affinity 
{cb = 10 ksT), intrinsic reaction attempt frequency {ko = 10^ s“^) and the diffusion 
constant {D =5pm^/s) is chosen such that the nucleation of domains and the unbinding 
of bonds are rare. Furthermore, we investigate the reaction and diffusion limited growth 
regimes for stiff (A = bksT/a~‘^) and soft receptors (A = 2kBT/a~‘^) mimicking bulky 
cell adhesion receptors and glycoprotein receptors, respectively. 
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Figure 10. Simulation results of the reaction limited radial growth for stiff receptors 
(left panels) and soft receptors (right panels). The first row (a, b) shows snapshots 
of the growing domain as a function of time whereas the second row (c, d) shows the 
number of bonds in the domain as a function of time. In the third (e, f) and the fourth 
row (g, h), we present bubble charts of the binding and unbinding rates depending on 
the number of neighboring bonds during the growth phase. The area of the bubbles 
represents the number of reactions with charts. 
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For ligand densities larger than the receptor density, we expect a qnadratic growth of 


the domain area 19,42,62 containing A'";, uniformly distributed bonds 

mt) = 


( 20 ) 


Here, po is the initial density of ligands and is the effective rate at the rim. 
This expression explicitly takes into account the two-dimensional nature of the growth 
process. 

The results of our Monte Carlo approach (blue full lines in Fig. 10 c and d) 
conhrm that the growth of the domain is, quadratic as expected. This is conhrmed 
by the excellent agreement of the data with £t by equation (20), shown in Fig. 10 


with dashed orange lines. The observed processes show that growth is faster for stiff 
{K'^=3.7 X 10"^ s“^) than for soft (iF“=2.0 x 10^ s“^) receptors, presumably because of 
stronger correlations between bonds. Clear deviations from the quadratic behavior take 
place when the finite size effects start to play a role and the domain begins to cover the 
whole simulation box. 

In order to relate the rates extracted from the £t to the microscopic rates which 
were actually used to grow the domains, we construct bubble charts for binding and 


unbinding rates (Fig. 10 e-h), which are classihed by the number of neighbors. A hxed 
number of neighbors can be organized in several different configurations around the 
receptor of interest, which results in the multiple bubbles for each number of neighbors. 
In the bubble charts, the area of the bubble is associated with the occurrence of a 
particular rate in the simulation. 

Interestingly, the effective rate K'^ corresponds very well to the average rate 
recorded in the simulation. Actually, for the stiff bonds the average rate at the 
rim, obtained by averaging all rates forming with up to hve neighbors K°^, is 
A'°°=3.7 X 10^s“^, and for soft bonds K°^=2.5 x 10^s“^. Rates for the formation of 
bonds with three to hve bonds in the neighborhood are most commonly observed (largest 
bubbles), which is consistent with the formation of new bonds at the edge of the domain. 
The rates for the formation of bonds with six or more neighbors are considered to be the 
results of events from rebinding within the domain, in agreement with the large number 


of dissociation events with seven and eight neighbors (Fig. 10 g, h) 


The analysis of microscopic rates in the bubble plots shows that the binding rates 
have a tendency to increase up to hve neighbors. This happens because the formation of 
additional bonds, in principle, reduces the distance between the receptor and the ligand 
at the position of the binding site. The rates for forming the bond with 3-5 neighbors 
are signihcantly larger for stih receptors, which is the source of the diherence in the 
speed of the overall growth process of the domain. The reason for this diherence is 
that for stih receptors, the membrane approaches closer to the substrate than for soft 
receptors, which themselves deform while forming a bond, leaving the membrane at a 
larger height. Rates for forming a bond with 6-8 neighbors decrease with increasing the 
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number of adjacent bonds. This effect is more significant for stiff receptors, because the 
fluctuations in the membrane are suppressed to a larger extent, and while the distance 
from the receptor is relatively small, stronger thermal membrane excitation is necessary 
to bring the ligand into the reaction zone of the receptor. The unbinding rates occur less 
frequently. The most common unbinding rate is the one with 8 adjacent bonds, which 
is clearly associated with unbinding within the domain. The unbinding rates decrease 
exponentially as a function of the number of neighbors, for both stiff and soft receptors, 
showing the stabilization effects that binding in the surrounding has on the respective 
bonds. 


6.3. Diffusion limited growth 


For ligand densities lower than the receptor density, the growth is diffusion limited (Fig. 
11 ), and depends only implicitly on the effective reaction rates through the density of 
ligands and bonds at the edge of the domain. In other words, the growth explicitly 
depends only on the diffusion constant, and the area of the growing domain A{f) is 


given by 19,21 


A{t) = Aira^Dt. (21) 

In this equation, a is a dimensionless speed factor, which is, in two dimensions, 
determined from the implicit equation 


PO — Pe 
Pb 


= exp Ei (a^). 


( 22 ) 


Here, pe is the density of ligands at the edge of the domain and Ei (x) is the so-called 
exponential integral 59 . This relation is obtained from the binder conservation at the 


rim of the domain and the respective solution of the diffusion equation (see Shenoy and 
Freund for details). 


We numerically solve equation (22) using the densities of bonds and ligands at the 


edge of the domain evaluated from the radially averaged density prohles. We obtain 
a = 0.34 for stiff receptors and a = 0.27 for soft receptors. The difference in the speed 
factors emerges from the somewhat larger density of bonds at the rim of a domain 
with soft receptors. Using these speed factors, we can calculate the expected diffusion 
constant from the linear £t (orange dashed lines in Fig. 0- Specihcally, we obtain 
a diffusion constant of 4.8pm^/s for the large bond stiffness, and 5.2pm^/s for the 


low bond stiffness (right column of figure 10) which is is consistent with the diffusion 
constant in the simulation (5.0pm^/s). 


6 . 4 . Remarks in the experimental context 

The obtained results from simulations of the growth of ligand-receptor domains show 
that it is, in principle, possible to relate macroscopic measurements with the underlying 
microscopic processes. From diffusion limited processes we can extract the diffusion 
constant with excellent accuracy, which is also accessible from experimental data. 
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Figure 11. Simulation results of radial growth in the diffusion limited case. Top: 
Snapshots of the growing domains. Bottom: Growth curves with linear fit indicating 
diffusion limited growth. Parameters like in Fig. [^except for initial ligand density 
(only 0.4 X 10® diffusing ligands). 


However, as noted before 42 , issues may arise if the crossover to the saturation of 


the growth curve due to the hnite size of the vesicle or cell occurs relatively quickly 
and the vesicle runs out of free binders. Furthermore, it is possible to relate the mean 
reaction rate to the microscopic events. 

7. Conclusions 


We presented two different approaches for simulating protein-mediated adhesion between 
membranes. The hrst simulation scheme considers the deformation and the fluctuations 
of the membrane explicitly, by evolving the membrane prohle with the help of a Langevin 
equation. The latter was derived from the Helfrich Hamiltonian and included the 
hydrodynamic interaction between membrane and surrounding fluid. The binding and 
unbinding of ligands and receptors is modeled by Dembo’s rates that are in detailed 
balance with the instantaneous shape of the membrane. Simpler variants of this scheme 
have been used successfully in earlier studies to describe thermal equilibrium and 
reaction limited dynamics 32 . However, this scheme fails to describe the dynamics 


on longer length scales as well as diffusion limited processes. The problem arises from 
the fact that time step is as short as 10“® s to correctly recover the membrane thermal 
excitations. Furthermore, the calculation of the membrane prohle requires the use of 
Fast Fourier Transformations which scale the simulation time with A^log(iV), where N 
is the number of considered membrane segments. As a result, only membrane patches 
of about pro? carrying about 1000 proteins can be simulated for about 0.1 seconds. 
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We overcome these constraints by constrncting an effective Monte Carlo scheme. 
In this scheme, we coarse-grain the adhesion dynamics by integrating the effects of the 
membrane into a set of effective reaction rates for ligand-receptor (nn)binding. These 
rates are derived by averaging Dembo’s rates over the membrane height flnctnations, 
which we do semi-analytically for an arbitrary bond conhgnration. This allows ns to 
circnmvent the explicit treatment of the membrane, and use a much larger time step 
in the simulation. Consequently, cell-sized objects (10"^/xm^) carrying 10® proteins can 
be simulated for several tens of seconds with the resolution of 10 nm and 10“® s. In 
this scheme, the simulation time scales linearly with the number of binders and the 
simulation time is thus reduced by a factor of about 10® for the parameters used in 
this study compared to the Langevin approach. As shown by an in-depth analysis of 
the correspondence between the Langevin and Monte Carlos simulations, this increased 
efficiency is achieved basically without loss of accuracy from the nucleation of adhesion 
domains and the early stages of growth to the asymptotic growth behavior and the 
saturation to an appropriate equilibrium state. 

This outstanding performance allows a successful study of completely realistic 
scaffolding processes. As an example, we performed an analysis for radially growing 
domains, which is one of the most common scenarios for the development of adhesions. 
We demonstrate that the measurables that can be extracted from the macroscopic 
development of the domain can be related to underlying microscopic stochastic 
processes, namely the protein diffusion and the binding kinetics. 

The simulations presented herein set a foundation for an in-depth analysis of protein 
transport and complexation dynamics in membranes, which is key to the understanding 
of the formation of functional microdomains and rafts. Furthermore, processes which 
present slow convergence or require correlations and signaling on the level of the entire 
cell are within the reach of accurate modeling. Now that the adhesion on the level of the 
membrane can be studied in great detail, the challenge becomes to couple the membrane 
to other cell structures and processes, which is a direction for future development. 
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Appendix A. Calculation of the membrane height distribution 

The membrane height distribution depends on the bond conhgnration of the membrane 
as well on the position of the binding site as can be seen in following equation. By 
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definition 

p{h{r)) = j P[/i'(r)]p[/i'(r)]5(/i'(r) - h{r)), (A.l) 

where we have on the left side the probability distribntion of the height p{h{r)) at 
the binding site r, whereas on the right side p[h'{r)] is the probability for having a 
membrane profile h'{r). This probability depends on the bond configuration (i.e. the 
positions ri of the ligand-receptor bonds). For simplicity, we set fd = {kBT)~^ = 1. 
To evaluate the above integral, the Boltzmann weight for p[h'{r)] determined by the 
Helfrich-Hamiltonian ([^ is plugged in and the Dirac function is written in the Fourier 
representation, which gives 


p{h{r)) (x du V[h{q)] exp 


Nb 


o 


i=l 


2A 


||h(q)|p + +iiy(h'(r) - h(r)) 


(A.2) 


We now apply successively N/, Hubbard-Stratonovich transformations, one for each bond 
term in the sum over i. This produces Nf, Gaussian integrals over auxiliary fields (pi. 
Furthermore, we write h{r) in the Fourier representation. As a result, we get 


p{h{r)) oc y dz/ y V[h{q)] j d(j)j 

X expj - ^ ^ ^ ilop)j - * ^ ^(q) exp (iq ■ rj) 


X] ll^(q)ir + 7 ) + ^ h(q) exp (iq ■ r) - iz//i(r) 


2 A 


Integration over h(q) is Gaussian integral and leads to 
p{h{r)) oc 


/ ( n / ( ~ XI ^ + X 

Vi / ^ i i 

exp (iq ■ rj)0j — z/exp (iq ■ r)^ 


Kq'^ + 7 


X 


exp (-iq ■ rk)0fc - v exp (-iq ■ r) 


(A.3) 


(A.4) 


In the following step, the terms within the curly brackets of equation (A.4) are 


reorganized, and the sums over q are converted to integrals that can be evaluated 
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independently leading to Kelvin functions (37). After some sorting, we obtain 

'5jk kei (go|rj - Fkl)' 


p{h{r)) (X I du in / 


A 


27rydf7 


=M, 


jk 


(A.5) 


j 


exp 


1 Z /2 


2 Sydfy 


— z/ z 


i/z(r) + ^( 


kei (go|rj - r |) 
' 27rydf7 


where 




(A.6) 


is the inverse of the lateral correlation length. Performing the Gaussian integrals in u 
gives after some algebra 


1 


p{h{r)) oc JJ / d0j exp - - 0* M. 


jk 


ijk 


2 kei (go|r - rjO^ei (go|r - fk) 




=Gjk{r) 


exp 


(^o + h{r 




TT 


(A.7) 


Since the remaining integrals are again Gaussian, one finally gets 

,4kei(go |r - rjl) 


p{h{r)) ocexp 


h(r)^/oGij(r) 


TT 




1 

2 


ij 


16kei(go k - ri|) Gij(r) ^ kei (go |r-rj|) 


TT^ 


h^(r) 


(A. 8 ) 


As the probability distribution (A. 8 ) is itself a Gaussian distribution, again, the average 
height can be calculated by completing the square in the exponent. Gonsequently, one 
obtains 


(h(r)) = h{v) = 


loGij (r) ^ kei (go I r - Fi I) 


Syifq + Ep kei (go |r - ri|)Gij(r)-i kei (go |r-rj|) 


(A.9) 


The fluctuations are simply given by 

(h^(r)) - {h{r)f = a(r) = 

— 16 kei (go |r - ri|) Gij(r)-^ kei (go |r-rj|) 


(A.IO) 
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